## pairwise comparisons distractor photo race correlates with photos observed
## only target photos used for analysis
##
## August 9, 2019
##
## Kevin Quinn
## University of Michigan
##

library(MCMCpack)
set.seed(592005)

mydata <- read.csv("../ScaleRaceSpring2017clean.csv")


## subset data
mydata <- mydata[mydata$condition %in% c("P1W", "P2B"),]

## keep responses and photo IDs but nothing else
target.inds <- c(6, 7, 10, 12, 15, 17, 19, 20, 22, 23)
keep.vars <- c(paste("X", target.inds, "_Q58", sep=""),
               paste("X", target.inds, "a.id", sep=""),
               paste("X", target.inds, "b.id", sep=""))

mydata.sub <- mydata[, keep.vars]



## reshape into long format
mydata.sub.long <- reshape(mydata.sub, direction="long",
                           varying=list(c(1:10), c(11:20), c(21:30)),
                           v.names=c("Y", "photoID.a", "photoID.b"),
                           ids=rownames(mydata.sub))

## keep only obs with non-missing values for MM
mydata.sub.long <- na.omit(mydata.sub.long)
mydata.sub.long <- mydata.sub.long[mydata.sub.long$Y != "",]


## convert to character variables
mydata.sub.long$Y <- as.character(mydata.sub.long$Y)
mydata.sub.long$photoID.a <- as.character(mydata.sub.long$photoID.a)
mydata.sub.long$photoID.b <- as.character(mydata.sub.long$photoID.b)

## recode Y to match what MCMCpaircompare expects
for (i in 1:nrow(mydata.sub.long)){
    if (mydata.sub.long$Y[i] == "the person in photo 1"){
        mydata.sub.long$Y[i] <- mydata.sub.long$photoID.a[i]
    }
    if (mydata.sub.long$Y[i] == "the person in photo 2"){
        mydata.sub.long$Y[i] <- mydata.sub.long$photoID.b[i]
    }    
}

mydata.sub.long <- mydata.sub.long[, c(5, 3, 4, 2)]

mydata.sub.long <- mydata.sub.long[sample(1:nrow(mydata.sub.long),
                                          size=1000, replace=FALSE),]

cat("\n\nN =", nrow(mydata.sub.long), "\n\n") 


## starting values and set up
pnames <- sort(unique(c(mydata.sub.long$photoID.a, mydata.sub.long$photoID.b)))

raw.counts <- table(mydata.sub.long$Y)
zero.names <- pnames[!(pnames %in% names(raw.counts))]
if (length(zero.names > 0)){
  old.names <- names(raw.counts)
  raw.counts <- c(rep(0, length(zero.names)), raw.counts)
  names(raw.counts) <- c(zero.names, old.names)
}
raw.counts <- raw.counts[pnames]
raw.ranks <- rank(raw.counts)
theta.start <- 2*(raw.ranks / length(pnames) - 0.5)



b2.pairwise.out <- MCMCpaircompare(mydata.sub.long,
                                   theta.constraints=list(mw045.J1="-",
                                                          mb037.J1="+"),
                                   theta.start=theta.start,
                                   alpha.fixed=TRUE,
                                   burnin=50000, mcmc=1000000, thin=50,
                                   verbose=10000, seed=24159
                         )


save(b2.pairwise.out, file="MM3.2.b2.pairwise.out.Rda")
